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Abstract. Given a set of basic areas, the territory design problem asks to create a predefined number of territories, 
each containing at least one basic area, such that an objective function is optimized. Desired properties of territories 
often include a reasonable balance, compact form, contiguity and small average journey times which are usually 
encoded in the objective function or formulated as constraints. We address the territory design problem by devel¬ 
oping graph theoretic models that also consider the underlying road network. The derived graph models enable us 
to tackle the territory design problem by modifying graph partitioning algorithms and mixed integer programming 
formulations so that the objective of the planning problem is taken into account. We test and compare the algorithms 
on several real world instances. 
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1 Practical Motivation 

The starting point for this research was a real territory design problem at the PTV AG; design 
and plan service territories for companies that send agents to customers on a regular basis. In this 
case, the number of territories to be planned is predetermined by the number of service agents. 
Required properties of territories are manifold. For example, the obtained territories should be 
reasonably balanced with respect to an activity index of the customers, they should have “good” 
accessibility (non-accessible areas have to be accounted for), have to be contiguous and compact; 
and most importantly the expected journey time for visiting a typical set of customers within the 
territory should be small. Even small improvements in the objective save fuel and time that is 
needed to visit the customers. In turn, the costs of the service company are reduced and the usage 
of valuable resources is optimized. This optimization can also lead to the company serving even 
more customers using the same amount of employees or increasing the face time with its customers. 

More formally, given a set of basic areas such as zip code areas or customer sites, the territory 
design problem asks to create a predefined number of larger territories such that an objective func¬ 
tion encoding the desired properties is optimized and territory constraints are fulfilled. Apart from 
the introduced example, other important applications include political districting, electrical power 
restricting or planning territories for social facilities like hospitals or administrative units. 

We organize the paper as follows. We begin in Section|^by introducing basic concepts and by 
summarizing related work. The main parts of the paper are Section and Section The former 
introduces our rationale, graph models and the approach employing graph partitioning algorithms 
to perform territory design. The latter uses a mixed integer programming approach to tackle the 
problem. We report on carefully designed experiments performed on several real world instances 
in Section Experiments indicate that our algorithms compute high quality territories that have 
small average journey times in particular. Einally, we conclude with Section 


2 Preliminaries 


2.1 Basic concepts 

Consider a set of basic areas B = {6i,..., with an activity index a : B ]R>o. We use the 
short notation a* for a{bi) and extend the activity index to sets, i.e. , for a set B' C B, its activity 
a(-B') := ^ territory T is defined as a subset of the set of basic areas. The activity 

index of a territory T is defined as a(T) = The territory design problem is to assign 

each basic area to one of the k territories, where /c > 1 is the number of desired territories given 
in advance. A solution to the territory design problem is a set of territories T = {Ti ,..., T*,} such 
that for each pair of territories, Tj fl = 0 and Ti U ■ ■ ■ U = 5. A balance constraint on the 
activity index demands that a(T) < (l+e)[^^] VT G T for some imbalance tolerance paxametev 
e > 0. A solution is called feasible if it obeys the balance constraint and other constraints placed 
on the contiguity and compactness of territories. The objective is to minimize the sum of pair-wise 
travel distances within the territories. Note, for small values of e a feasible solution may not exist. 

Consider an undirected graph G = {V,E,c,u) with edge weights u : E ^ M>o and node 
weights c : V —>■ M>o. Again, we extend c and u to sets, i.e. , for a set V' <fV its node weight 
c{V') := Yhv&v ^('^) ^ E' (Z E its edge weight u:{E') := Yhe&E' The set r{v) : = 

{u : {v,u} G E} denotes the neighbors of v. A maximal connected component is a maximal 
subgraph of G in which every pair of nodes are connected by a path. The graph partitioning 
problem is looking for k blocks of nodes Vi,... ,14 that partition V, i.e., 14 U ■ ■ ■ U 14 = V and 
ViHVj = 0 for z 7 ^ j. Here, a balancing constraint demands that Vz G {1,..., A;} : c(14) < Tmax : = 
(1 + e) \c{V)/k~\ for some imbalance parameter e > 0. Often the objective of graph partitioning 
problems is to minimize the total cut where Eij := {{u,v} e E : u E Vi, v e Vj}. 

However, during the course of the paper we will modify this objective. 

2.2 Related Work 

There has been a huge amount of research on territory design so that we refer the reader to the 
surveys mm for most of the material. For material related to graph partitioning, we refer the 
reader to the survey [|2]. Here, we focus on issues closely related to our main contributions. 

The basic formulation of territory design problems is due to Hess and Samuels BH who used 
relaxations of integer linear programming models to tackle the problem. Ronen [fTOll presents a 
mixed integer formulation to minimize the total driving distance of salesmen. However, it is as¬ 
sumed that each trip to a customer requires a separate trip which includes that the salesmen returns 
home, so that the average journey time is not measured. 

Using the same assumption, Zoltner and Sinha lITSl develop a linear optimization model incor¬ 
porating road networks for the distance computation to minimize the travel time. Additionally, the 
networks are used to represent connectivity between the territories which helped to compute con¬ 
tiguous territories with improved accessibility. Note that in both cases the objective also implicitly 
encodes compactness of the territories. 

Forman and Yue [[5l use a genetic algorithm to compute congressional territories. The basic idea 
is based on an encoding and on genetic operators that where originally used to solve the Traveling 
Salesman Problem. Desired properties of the territories are integrated into the fitness function of 
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the individuals. In this paper, we also use an evolutionary algorithm whieh instead eneodes the 
individuals as partitions of a graph model. 

Hess et al. ||71 deseribe a mixed integer program based on the loeation-alloeation approaeh. In 
eaeh iteration basie areas are assigned to territory eenters and afterwards the eenters are updated. 
Compaet territories are aehieved by minimizing the sum of the distanees of the basie areas to the 
territory eenters. Sinee we also implement and extend this approaeh we go into more detail later. 

Reeently, Butseh et al. [|3l used a reeursive geometrie bipartitioning approaeh to assign basie 
areas to territories. The approaeh reeursively eomputes bipartitions of the basie areas using their 
eoordinates. However, this approaeh may eompute longish, non-eompaet territories and has prob¬ 
lems ineorporating geographie obstaeles sueh as rivers or mountains. We eompare our algorithm 
against this approaeh in Seetion]^ 

2.3 KaHIP 

Within this work, we use the open souree multilevel graph partitioning framework KaHIP [IIIII3II 
(Karlsruhe High Quality Partitioning^^ More preeisely, we modify the distributed evolutionary 
algorithm KaFFPaE eontained therein to ereate partitions of our graph models whieh in turn yield 
solutions to the territory design problem. Henee, we shortly outline the main eomponents of KaHIP. 

Besides the evolutionary algorithm, KaHIP implements many different algorithms, for exam¬ 
ple flow-based methods and more-loealized loeal searehes within a multilevel framework ealled 
KaFFPa, as well as several eoarse-grained parallel and sequential meta-heuristies. The algorithms 
in KaHIP have been able to improve the best known partitioning results in the Walshaw Beneh- 
mark [[T4l for many inputs using a short amount of time to ereate the partitions. 

Evolutionary Algorithm Outline. We now roughly outline the general strueture of KaFFPaE sinee 
one of our algorithms employs a modified version of the evolutionary algorithm. KaEEPaE starts 
with a population of individuals (in our ease partitions of the graph) and evolves the population 
into different populations over several rounds. In eaeh round, the evolutionary algorithm uses a 
seleetion rule based on the fitness of the individuals of the population to seleet good individuals 
and eombines them to obtain improved 


Algorithm 1 General Stueture of the Evolutionary Graph Partitioning Algorithm KaEEPaE. 
ereate initial population P 
while stopping eriterion not fulfilled 

select parents pi, P 2 from P (using fitness funetion) 

combine pi with p 2 to ereate offspring o (using a graph partitioner) 

mutate offspring o 

evict individual in population using o 
return the fittest individual that oeeurred 


* available at http : //algo2 . iti . kit. edu/documents/kahip/ 
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offspring. More precisely, KaFFPaE uses a tournament selection rule to select individuals for com¬ 
bination. In KaFFPaE, the fitness function is set to the number of edges cut, however, we will 
modify the fitness function for the purpose of territory design. A combine operation then employs 
a graph partitioning framework to obtain an offspring having the good cuts of both input partitions. 
We refer the reader to IfT^ for more details. Our algorithm generates one offspring per generation. 
The general structure of the evolutionary algorithm is depicted in Algorithm 

The algorithm is parallelized by giving each processing element its own population so that 
combine and mutation operations can be performed independently. This is combined with a scal¬ 
able communication protocol to exchange high quality solutions between the processing elements 
over time. 

3 Territory Design by Graph Partitioning 

The territory design problem and the graph partitioning problem are closely related. Our approach 
to territory design using graph partitioning consists of two steps. At first we construct a graph 
that corresponds to the territory design problem. Then, a customized algorithm partitions the con¬ 
structed graph into k blocks. We develop a one-to-one mapping between basic areas and vertices 
in the constructed model. Hence, a partition of our model yields a solution of the territory design 
problem. Edges are defined by a neighboring relation. We begin this section by explaining how we 
construct the graph to be partitioned, then define the fitness function that is used in the evolutionary 
algorithm and outline how everything is put together. 

3.1 Constructing the Graph 

We now explain how we construct the graph that will be partitioned by the graph partitioning 
algorithm which in turn gives us a solution to the territory design problem. Every basic area bi with 
activity Oj corresponds to a node Vi in the graph with weight q := Oj. An edge between a node Vi 
and Vj exists if and only if the basic areas bi and bj are neighbors as described below. In general, 
we set the weights of the edges in our model to one. However, as we will see later, the evolutionary 
graph partitioning algorithm that is used to partition the model takes the real distances into account 
to compute the value of the objective function of a solution. 

Edges in our Model. The algorithm we use to compute the edges in our model is as follows. Note 
that we have multiple design goals for our model. First of all, two basic areas that are close should 
be connected by an edge in the model. Additionally, the model should not be too dense, e.g. if the 
maximum degree of the graph is bounded, the graph should be connected and edges that are too 
long should also not be contained in the model. Roughly speaking, we perform two iterations of 
Kruskal’s algorithm 0 on a complete graph where every node Vi corresponds to a basic area bi 
and the edge weights are the distances between the basic areas uj{vi, Vj) = d{bi, bj). Recall, that 
Kruskal’s algorithm scans the edges of the graph in increasing order of their weight to grow a forest 
and that it adds the edges joining two trees in the forest. 

After the first run of Kruskal’s algorithm, every edge that is in the minimum spanning tree 
(MST) computed by the algorithm is inserted into our model. For the second run of Kruskal’s 
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algorithm, we remove the just eomputed MST edges from the eomplete graph. Now, while the 
algorithm seans the edges in inereasing order, it adds the eurrent edge to our model if the maximum 
node degree in the eurrent state of the model does not exeeed a user defined parameter 7 > 1 and 
if the length of the edge is smaller or equal to /3 ■ cuavg- Here, /? is a user speeified faetor and Wavg 
denotes the average edge weight of MST edges of the first iteration of the algorithm. 

3.2 Fitness Function 

Reeall that we are looking for a partition of the just defined graph model into k bloeks of nodes 
P = {Vi,..., 14 }. Moreover, eaeh subgraph indueed by a bloek Vi should to be contiguous. In 
other words, if rimcciVi) is the number of maximal eonneeted eomponents in the subgraph indueed 
by 14 , we want rimcciVi) = 1 for eaeh bloek VJ. We define the number of eonneeted eomponents of a 
partition P as ncon{P) = ^mcc(Vi). To guide the evolutionary algorithm towards eontiguous 
bloeks, we use a penalty approaeh. More preeisely, we use the faetor (1 + a{ncon{P) — k)) with 
penalty parameter a > 0 in our objeetive funetion to ensure that non-eontiguous territories are 
penalized. We then set the objeetive/fitness funetion of the evolutionary algorithm that partitions 
our model to 

k 

fobj{P) = (1 + a{ncon{P) - k)) 

i=l 

where d{u, v) is the average time needed to traverse the shortest route in the original network 
between the basie areas eorresponding to u and v respeetively. 


3.3 Overall Algorithm 


To taekle the territory design problem, we use the distributed evolutionary graph partitioner KaFF- 
PaE to partition our model. We modify the fitness funetion of the algorithm to the objeetive funetion 


presented in Seetion 3.2 Note that the partitioning algorithms within KaFFPaE still optimize the 
number of out edges. The amount of allowed imbalanoe is a oonstraint of the partition problem in 
KaEEPaE. Due to the looal seareh algorithms in KaEEPaE no bloek of a partition will be empty. 
However, it may be possible that partitions oreated during the oourse of the evolutionary algo¬ 
rithm are not eontiguous. Henoe, whenever we oreate an individuum/partition, we try to make it 
eontiguous. This is done by grouping neighboring eonneeted eomponents of the partition so that 
eaeh bloek beeomes eontiguous. More preeisely, exeess eonneeted eomponents are assigned to the 
neighboring bloek with the least aetivity. Doing this may result in imbalaneed partitions so that 
a rebalaneing step is performed afterwards to ensure the balanee eonstraint. We eall the overall 
algorithm to taekle the territory design problem KaTeD (Karlsruhe Territory Design). 
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4 Territory Design by Mixed Integer Programming 

Besides the graph partitioning approaeh, we use a loeation-alloeation method for territory design. 
In this seetion, we first deseribe the loeation-alloeation method introdueed by Hess et al. iQ and 
then present the details of our modified and extended loeation-alloeation approaeh. 


4.1 The location-allocation approach hy Hess et al. 

The idea of solving territory design problems by means of a loeation-alloeation approaeh goes baek 
to Hess et al. Q. They apply the approaeh to the design of legislative distriets. In this applieation, 
territories are supposed to be compaet, contiguous, and balanced in terms of population. Due to the 
complexity of the problem, Hess et al. DTI decompose it into two subproblems which are solved in 
an iterative manner: 


- The location subproblem seeks to find (virtual) territory centers which are used to calculate the 
compactness measure. In the first iteration, the centers are a guess; in all subsequent iterations 
they are obtained by computing the center of gravity in each territory. 

- In the allocation subproblem the basic areas are assigned to territory centers. To this purpose, 
Hess et al. [Tl formulate an integer linear program. Let C denote the set of territory centers 
and dh-c- the distance between basic area bi E B and center cj E C. Furthermore, define the 
following decision variables: 




1 if basic area bi is assigned to center Cj 
0 otherwise 


Then, the model of Hess et al. Q can be stated as follows: 


min Y dlc.aiXb^cj 

Cj^C 


(1) 

Y/ ^ ^ 

CjGC 

bi E B 

(2) 


Cj E C 

(3) 

Xbicj e {0,1} 

bi ^ B ^ Cj ^ C 

(4) 


The objective function ([T]) optimizes compactness which is measured as the sum of the squared 
distances between basic areas and associated territory centers, weighted by the basic areas’ 
activity index. Constraints Q in combination with the integrality conditions Q ensure that 
each basic area is assigned to exactly one territory center. Balance is achieved by the constraints 
in Q which guarantee that the activity index of a territory deviates by at most e ■ 100 percent 
from the mean activity index. Instead of solving the integer program, Hess et al. Q set e to 
zero, solve the linear programming relaxation and then resolve all fractional assignments. 

Hess et al. [[Til perform location and allocation alternately until the solution converges. The algo¬ 
rithm repeats if multiple initial guesses for the territory centers are available. 
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4.2 Modifications and Extensions 

We now outline our modifications and extensions of the approach by Hess et al. [|71|. Roughly 
speaking, we use a well-known procedure to determine good initial centers, modify the balance 
constraint, solve the integer program directly and apply a multi-start procedure. We call our version 
of that approach KaLocAlloc (Karlsruhe Location Allocation). 

Location Step. The quality of the solutions obtained by the location-allocation approach strongly 
depends on the selection of the initial centers, i.e. , the centers used in the first iteration of the 
algorithm. A good initial set of centers should be well-distributed across the region under study. 
To this end, we adopt the initialization procedure from the k-means-i-i- algorithm [HI : Among all 
basic areas we pick the first center uniformly at random. Among all remaining basic areas the next 
center is picked with a probability which is proportional to the basic area’s squared distance to the 
nearest center already chosen. This is repeated until k centers have been selected. 

Allocation Step. For the allocation of basic areas to centers we use model ([T]) - Q with one 
modification. In order to limit only the maximum activity index of the territories, we drop the 
lower bound in constraints 


diXbiCj < (1 + 'i Cj eC ( 5 ) 

bi£B 

Although we drop the lower bound, solutions cannot contain empty territories. This is because 
the centers are picked among the basic areas and, therefore, at least the basic areas corresponding 
to centers are assigned to the respective center due to a distance of 0. Another difference to Hess 
et al. Q is that we do not solve the linear programming relaxation of the model, but the integer 
program. 

Multi-start Procedure. As already mentioned, the selection of good initial centers is important 
to achieve solutions of high quality. Therefore, we apply a multi-start procedure in which the 
problem is solved multiple times with different initial centers at each start. The centers at each 
start are selected randomly according to the k-means-i-i- scheme described above. After a user- 
defined number of starts the approach returns the best solution across all starts. Since in our case 
the objective is to minimize the sum of pairwise travel times between all basic areas of the same 
territory, the multi-start algorithm returns the best solution according to this criterion. 
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5 Experimental Evaluation 


Methodology. All experiments were performed on real-world data provided by PTV Group. The 
test data eomprises 15 instanees whose sizes range from approximately 300 to 5,000 basie areas. 
The number of territories to be planned is given in the test data and varies from 3 to 46. Further¬ 
more, road distanees and travel times for the shortest path between eaeh pair of basie areas were 
available and have been used for eompaetness evaluation. We also eompare our algorithms against 
the reeursive partitioning approaeh by Butseh et al. [|3l (BKNS) whieh has been provided by the 
authors. However, we perform only one repetition sinee the algorithm is deterministie. 

We evaluated our approaehes on all 15 instanees with a time limit of 300 seeonds. The im- 
balanee toleranee parameter e was set to 0.05 in all experiments. From the praetieal experienee of 
PTV Group the ehosen values for the time limit and for the imbalanee toleranee are very aeeeptable 
values for human planners. 

Parametrization. Parameters whieh are speeifie to territory design by graph partitioning were set 
as follows: penalty parameter to inerease eonneetedness a = 0 . 1 ; edge length faetor in our model 
= 5; node degree bound in our model 7 = 20. We repeated both approaehes 40 times using 
different random seeds for initialization and the average was taken. Parameters whieh are speeifie 
to the loeation-alloeation approaeh were set as follows: For the alloeation step, the relative MIP 
optimality gap was set to 0.001. The maximum time spent on the alloeation step was limited to 15 
seeonds, and the number of multi-starts was unrestrieted. 

System. All experiments have been done on an Intel Xeon CPU E3-1245 at 3,4GHz having 16 GB 
RAM and Mierosoft Windows 7. The loeation-alloeation approaeh was implemented in Java. We 
used Gurobi 5.6 to solve the integer linear program within the loeation-alloeation approaeh. The 
number of threads was set to 4 for KaTeD and KaLoeAlloe to ensure eomparability. 

5.1 Computational Results 

We shortly summarize the main results and present detailed per instanee results in Table First 
of all, the balanee eonstraint is satisfied in all eases for all approaehes. In 8 of 15 instanees, the 
average solution quality of KaTeD outperforms KaLoeAlloe. Considering the other 7 instanees, 
the value of the objeetive of the territories eomputed by KaLoeAlloe is 1 .6 pereent smaller than the 
ones eomputed by the evolutionary approaeh. The reeursive partitioning approaeh BKNS always 
yields worse results than both of our approaehes, but is faster sinee we eould perform only one 
repetition of the algorithm due to the faet that the algorithm is deterministie. On the largest instanee, 
G02, the BKNS algorithm needed 203 seeonds to eompute a result. On average, BKNS yields 
3.7 pereent larger objeetives than KaTeD and 4.3 pereent larger objeetives than KaLoeAlloe. The 
largest improvement over BKNS is obtained on instanee BL16 (Thiiringen) and amounts to 10.3 
pereent eompared to the result of KaTeD whieh eomputes the best objeetive on that instanee. 
Ligures to eompare the visual result eomputed by the different algorithms on two exemplary 
instanees, Thiiringen and Baden-Wiirtemberg in Germany. 
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6 Conclusion and Future Work 


In this paper we addressed the territory design problem by developing graph theoretie models that 
also eonsider the underlying road network. The derived graph models enabled us to taekle the 
territory design problem by redueing it to a graph partitioning problem. The resulting graph parti¬ 
tioning problem is then solved by using a modified evolutionary graph partitioning algorithm takes 
the objeetive funetion of the territory design problem into aeeount. On the other hand we extended 
an existing mixed integer programming formulation. We tested and eompared the algorithms on 
several real world instanees. 

Important future work ineludes the integration of the objeetive funetion of the territory design 
problem direetly into a multi-level graph partitioning algorithm. In partieular, it would be interest¬ 
ing to define loeal seareh algorithms for our objeetive. On the other hand, it would be good to have 
a graph partitioning algorithm that ean ensure eonneetedness of bloeks. 

Acknowledgements. We would like to thank Alexander Butseh for providing us with an implemen¬ 
tation of the reeursive partitioning approaeh of Butseh et al. [|3l . 
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A Tables and Pictures 


Table 1. Sum of Pairwise Travel Times - Average 


Instance 

#Basic Areas 

#Territories 

KaTeD avg 

KaLocAlloe avg 

BKNS avg 

BLl 

520 

6 

56 080 176 

54157925 

56943 136 

BL3 

1265 

13 

143 300619 

144435173 

146 538 257 

BL5 

2440 

20 

239 489 621 

238965327 

250332974 

BL6 

514 

6 

109625 833 

110561 195 

118 085 932 

BL7 

1038 

9 

48151371 

48 456 559 

50664649 

BL8 

1252 

11 

159771590 

160239 521 

172659 265 

BL9 

2019 

18 

260762 258 

258027628 

271895 453 

BLll 

427 

3 

28 990949 

28796550 

29745 445 

BL13 

289 

3 

50774924 

50326699 

50301370 

BL14 

493 

5 

55 345 856 

55 723 702 

55 844700 

BL15 

284 

3 

41953 610 

40997882 

41 426 146 

BL16 

428 

5 

44396659 

44564384 

48 977 328 

GOl 

4472 

41 

469 430937 

459630225 

479442432 

G02 

4971 

45 

597 515 028 

591514229 

632623 485 

G03 

2 241 

20 

333 612833 

329 314024 

343 170613 


Table 2. Sum of Pairwise Travel Times - Minimum and Maximum 


Instanee 

KaTeD min 

KaLocAlloe min 

KaTeD max 

KaLoeAlloe max 

BKNS min/max 

BLl 

54918613 

53663743 

56 276 617 

54819237 

56 

943 136 

BL3 

142464954 

143 406994 

143688813 

145 504 176 

146 

538 257 

BL5 

235 752526 

233110996 

244426840 

245 775 302 

250 

332974 

BL6 

109140977 

109 548 283 

109893677 

111336978 

118 

085 932 

BL7 

47118442 

48 293 411 

48 591642 

48569098 

50 

664 649 

BL8 

158613738 

158257253 

160899255 

161325 525 

172 

659 265 

BL9 

257 234587 

253678835 

262697 498 

261954880 

271 

895 453 

BLll 

28 973 547 

28796550 

28 992672 

28796550 

29 

745 445 

BL13 

50774924 

50217885 

50774924 

50446097 

50 

301 370 

BL14 

54102563 

55 456757 
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Fig. 1. Visualization of the KaTeD result of instance BL16 (Thiiringen, Germany). 



Fig. 2. Visualization of the KaLocAlloc result of instance BL16 (Thiiringen, Germany). 
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Fig. 3. Visualization of the BKNS result of instance BL16 (Thiiringen, Germany). 



Fig. 4. Visualization of the KaTeD result of instance BL8 (Baden-Wuerttemberg, Germany) 
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